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Abstract 

Army  installations  contain  a  high  density  of  mission-critical  buildings  that 
require  constant  protection  from  accidental  or  deliberate  contamination  of 
water  distribution  systems.  Current  simulations  of  contaminant  fate  and 
transport  in  pipe  systems  do  not  accurately  portray  reality.  The  simula¬ 
tions  assume  pure  hydraulic  transport  of  contaminants  and  do  not  account 
for  sorption  of  the  contaminant  on  pipe  walls.  Additionally,  subsequent 
reactions  such  as  hydrolysis  are  not  considered.  These  omissions  reduce 
predictability  of  a  contaminant’s  progression  and  effect  in  the  distribution 
system.  In  addition,  inadequate  understanding  of  contaminant  fate  and 
transport  may  result  in  an  unacceptable  risk  to  mission  readiness.  Howev¬ 
er,  performing  laboratory  tests  for  every  known  and  emerging  chemical  or 
biological  contaminant  to  obtain  uptake  and  reaction  parameters  is  not 
feasible  with  regard  to  time  or  cost  investments.  This  report  documents 
advances  in  molecular  modeling  predictions  for  the  transport  of  contami¬ 
nants  using:  the  Nanoscale  Molecular  Dynamics  program,  the  influence  of 
bacterial  biofilms  on  spore  viability  within  a  chlorinated  water  distribution 
system,  the  computational  chemistry  predictions  of  the  rate  of  hydrolysis 
of  a  specific  contaminant,  and  the  inclusion  of  results  into  the  predictive 
software,  EPANet.  This  work  opens  the  way  for  better  vulnerability  as¬ 
sessments,  protection,  real-time  response,  remediation,  and  planning. 


DISCLAIMER:  The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication,  or  promotional  purposes. 
Citation  of  trade  names  does  not  constitute  an  official  endorsement  or  approval  of  the  use  of  such  commercial  products. 
All  product  names  and  trademarks  cited  are  the  property  of  their  respective  owners.  The  findings  of  this  report  are  not  to 
be  construed  as  an  official  Department  of  the  Army  position  unless  so  designated  by  other  authorized  documents. 
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1  Introduction 

1.1  Background 

Army  installations  have  a  high  density  of  critical  facilities  that  are  difficult 
to  abandon  in  place,  especially  during  force  projection  and  deployment. 

Examples  of  such  facilities  include:  Combatant  Command  (COCOM) 
Headquarters,  Department  of  Defense  (DoD)  Command  and  Control  Cen¬ 
ters  (C2),  and  all  medical  facilities.  In  the  event  of  chemical-biological 
(CB)  contamination  of  the  water  distribution  system,  such  mission-critical 
facilities  and  their  supporting  installations  would  need  protection.  Accu¬ 
rate  understanding  of  the  fate  and  transport  of  contaminants  is  vital  to  en¬ 
suring  adequate  protection  from  intentional  or  accidental  contamination 
events. 

Previous  U.S.  Army  Engineer  Research  and  Development  Center  (ERDC) 
research1  had  formulated  a  physicochemical  model  that  described  how 
contaminants  introduced  to  a  potable  water  distribution  system  were  ca¬ 
pable  of  uptake  (sorption)  on  pipe  materials.  This  physicochemical  model 
was  validated  using  laboratory  measurements: 

^^-  =  /?(l-A'/r)  Equation  1 

Q 


where: 

P(t)  =  concentration  of  contaminant  on  pipe  material 
Co  =  initial  concentration  within  the  solution 
P  =  fraction  deposited  on  the  wall  at  equilibrium 
T  =  characteristic  time  constant 

The  two  parameters,  P  and  T,  are  measured  experimentally  for  each  com¬ 
bination  of  pipe  material  and  contaminant.  Below  the  saturation  point  of 
the  contaminant  in  water,  the  parameter  P  is  insensitive  to  concentration. 
This  finding  shows  that  the  process  is  dominated  by  thermodynamic  parti¬ 
tion  of  the  contaminant  between  water  and  pipe  material,  and  that  P  can 


1  Direct-funded  6.2  research  unit,  “New  Physicochemical  Models  of  Fate  and  Transport  of  Contaminants 
in  Potable  Water  Distribution  Systems." 
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be  measured  in  percentage  uptake.  This  is  a  major  update  to  physicochem¬ 
ical  modeling  of  fate  and  transport  of  contaminants.  Leaving  these  two  pa¬ 
rameters  ((3  and  T)  as  unknowns  that  must  be  measured  on  a  case-by-case 
basis  for  each  pipe-wall-material/contaminant  pair  is  not  a  totally  satisfy¬ 
ing  result.  In  that  scenario,  conducting  costly  lab  experiments  on  all  likely 
contaminants  that  could  interact  with  all  common  pipe  materials  ahead  of 
an  event  is  too  expensive,  and  conducting  the  same  types  of  experiments 
during  or  after  an  event  is  too  slow  to  offer  satisfactory  remediation. 

This  pitfall  was  known  at  the  time  of  ERDC’s  previous  research  so  that 
contaminant  data  was  also  compared  to  simulant  data  in  the  hope  that  up¬ 
take  could  be  predicted  ab  initio  by  using  other  commonly  known  proper¬ 
ties  of  a  prospective  contaminant.  The  results  of  this  line  of  research 
showed  that  chemical  contaminant  uptake  is  a  function  of  the  parameters 
given  below. 

•  Molecular  weight:  Large  molecules  would  have  difficulty  diffusing 
into  the  molecular  matrix  of  the  pipe  wall  material. 

•  Hydrophobicity:  Hydrophilic  materials  are  unlikely  to  have  any  af¬ 
finity  to  pipe  wall  material. 

•  Number  of  contact  points:  Increasing  the  number  of  contact  points 
between  the  contaminant  and  the  pipe  wall  leads  to  greater  uptake. 

At  the  outset,  the  first  two  parameters  were  known  to  researchers.  The  ex¬ 
perimental  results  obtained  during  this  project  guided  us  to  the  third  pa¬ 
rameter,  the  concept  of  “number  of  contact  points.”  By  obtaining  a  near 
match  for  both  hydrophobicity  and  molecular  weight,  a  high-quality  simu¬ 
lant  (simulant  P)  for  an  agent  under  study  (agent  V)  could  be  selected. 

Considering  the  prohibitive  cost  of  performing  laboratory  studies  on  the 
behavior  of  all  possible  contaminants  within  a  water  distribution  system, 
this  project  focused  on  the  use  of  computational  methods  for  the  ab  initio 
determination  of  rate  constants. 

This  report  demonstrates  the  inclusion  of  results  from  molecular  modeling 
for  the  transport  of  contaminants  by  using  the  NAnoscale  Molecular  Dy- 
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namics  program  (NAMD)1,  the  influence  of  bacterial  biofilms  on  spore  vi¬ 
ability  in  the  presence  in  chlorinated  water  distribution  systems,  and  com¬ 
putational  chemistry  predictions  of  the  rate  of  hydrolysis  of  a  specific  con¬ 
taminant  into  predictive  software  EPANet2.  The  inclusion  of  these  initial 
values  (Section  1.1)  opens  the  way  for  further  water-distribution  simula¬ 
tion  upgrades  for  better  protection,  remediation,  and  planning  ahead  of  a 
contamination  event. 

1.2  Objectives 

The  objective  of  this  work  is  to  implement  improvements  to  a  simulation 
tool  for  planning  courses  of  action  (COA),  passive  protection,  active  pro¬ 
tection,  and  decontamination.  This  is  done  by  using  experimentally  vali¬ 
dated  molecular  models  for  CB  contaminant  sorption/desorption  kinetics 
using  computational  chemistry  (CC)  and  NAMD,  resulting  in  a  pipe  seg¬ 
ment  (PS)  hydraulic  model  based  on  key  kinetic  parameters  and  variables. 

1.3  Approach 

The  approach  of  this  project  can  be  broken  down  into  three  tasks  that  were 
concurrently  executed: 

•  Task  l :  For  bulk  reactions  with  chlorinated  water  species,  the  equa¬ 
tions  built  into  EPANet  were  adequate;  however,  there  is  no  lookup 
mechanism  to  bring  in  the  correct  kinematic  parameters  into  the  simu¬ 
lation  based  on  the  contaminant  selected.  To  accomplish  this,  we  add¬ 
ed  routines  to  manage  a  library  of  such  parameters.  The  use  of  this  li¬ 
brary  should  be  transparent  to  the  user  who  merely  selects  the 
contaminant  of  interest. 

•  Task  2 :  For  pipe  wall  sorption  /  desorption,  the  equations  built  into 
EPANet  were  inadequate  and  we  added  them  to  the  software.  In  addi¬ 
tion,  these  equations  require  access  to  a  library  of  kinematic  parame¬ 
ters  that  is  quite  similar  to  that  developed  in  task  l. 


1  NAMD  is  a  parallel  molecular  dynamics  code  designed  for  high-performance  simulation  of  large 
biomolecular  systems.  NAMD  is  distributed  free  of  charge  with  source  code 
(www.ks.uiuc.edu/Research/namd/). 

2  EPANet  is  an  EPA-developed  (Water  Supply  and  Water  Resources  Division)  software  that  models  water 
distribution  piping  systems.  It  is  a  Windows-based  program  that  performs  extended-period  simulation 
of  the  hydraulic  and  water  quality  behavior  within  pressurized  pipe  networks 

http://www.epa  .gov/n  rm  rl/wswrd/dw/epa  net.  htm  I). 
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•  Task  3 :  The  kinematic  parameters  that  populate  the  libraries  were  de¬ 
rived  from  a  combination  of  molecular  models  (as  outlined  in  Section 
1.1  above)  including:  NAMD,  Cosmotherm,  and  Gaussian.  Also,  careful 
characterization  of  biofilms  was  important  from  the  standpoint  of  ask¬ 
ing  how  much  biofilms  contribute  to  sorption  and  desorption  and  if 
they  are  suitable  hosts  for  biological  species  of  interest. 

1.4  Scope 

The  project’s  scope  is  a  proof-of-concept  for  using  computational  chemis¬ 
try  to  describe  sorption-desorption  of  a  family  of  contaminants  grouped  by 
molecular  weight,  and  an  initial  modification  of  existing  hydraulic  model 
code  (EPANet  2.0)  to  account  for  sorption-desorption. 

1.5  Mode  of  technology  transfer 

The  technologies  developed  by  this  work  will  be  transferred  by  imbedding 
the  fate  and  transport  kinetic  parameters  obtained  into  EPANet  2  MSX, 
the  establishment  of  a  Memorandum  of  Understanding  (MOU)  with  the 
U.S.  Environmental  Protection  Agency  (USEPA)  to  continue  joint  research 
efforts  pertaining  to  water  security,  and  recommended  changes  to  DoD 
guidelines. 
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2  Experimental  Program 

2.1  NAMD  modeling 

Conceptually,  the  parameters  that  need  to  be  derived  from  NAMD  and 
then  plugged  in  to  EPANet  are  easily  summarized  from  the  standpoint  of 
statistical  mechanics.  We  start  by  quoting  the  results  required  from  the 
standpoint  of  thermodynamics  and  then  show  how  the  thermodynamic 
parameters  are  derived  from  NAMD  simulations. 

2.1.1  Thermodynamic  derivation  of  the  uptake  parameter  0  for  EPANet 

First,  we  show  that  the  uptake  parameter  (P)  is  directly  related  to  chemical 
potential.  In  the  range  of  concentration  above  detectable  limits  and  below 
saturation,  the  best  approximation  is  a  first-order  interaction  between  the 
contaminant  and  the  pipe  surface  that  relaxes  toward  the  equilibrium  con¬ 
stant  (3.  (One  way  of  viewing  this  interaction  from  a  conceptual  standpoint 
is  to  compare  it  to  the  widely  used  octanol-water  ratio.  Another  valid  way 
to  view  our  work  is  to  think  of  this  as  a  pipe-water  ratio.)  This  relationship 
is  reproduced  again  in  Equation  2  as  an  uptake  with  one  time  constant  and 
one  equilibrium  constant. 


Pit ) 


Equation  2 


where 

P  (t)  =  concentration  of  contaminant  on  pipe  material  as  of  a 
function  of  time 

P  =  fraction  deposited  on  the  wall  at  equilibrium 
e  =  Euler’s  constant 

Co  =  initial  concentration  within  the  solution 
T  =  characteristic  time  constant 


The  equilibrium  of  this  equation  is  characterized  by  letting  time  go  to  in¬ 
finity,  yielding: 


Equation  3 
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The  chemical  potential,  denoted  as  E  in  Figure  l,  is  the  amount  of  energy 
either  gained  or  lost  when  a  molecule  of  contaminant  moves  from  the  wa¬ 
ter  onto  the  pipe.  Chemical  potential  is,  in  essence,  the  binding  energy  of  a 
contaminant  to  the  pipe  surface. 


E, 


'Pipe 


^  ^Pipe"^AVater 


Ewater 


Figure  1.  Graphic  of  chemical  potential  -  the  binding  energy 
of  a  contaminant  to  the  pipe  surface. 

Hence,  the  equilibrium  state  of  the  system  is  characterized  by  Equation  4. 


Pit) 


C(t ) 


fi 


\-p 


Equation  4 


where 

P  (0  =  concentration  of  contaminant  on  pipe  material  as  a  function 
of  time 

C  (t)  =  concentration  of  contaminant  in  the  water  as  a  function  of 
time 

P  =  fraction  deposited  on  the  wall  at  equilibrium 

This  characterization  means  that  the  asymptotic  uptake  in  the  pipe  (/?)  is 
directly  related  to  chemical  potential.  If  the  contaminant  ends  up  almost 
entirely  on  the  pipe  or  almost  entirely  in  the  water  implies  that  |E|  >  >kT. 

Because  NAMD  is  computationally  efficient  enough  to  measure  the  bind¬ 
ing  energy  over  a  large  number  of  molecules  and  binding  events  with  the 
simulated  pipe  wall  surface,  we  arrived  at  a  valid  estimation  of  the  chemi¬ 
cal  potential  by  calculating  the  fraction  of  molecule  bound  to  the  wall  and 
claiming  that  this  calculation  corresponds  to  the  asymptotic  uptake. 

2.1.2  Derivation  of  the  uptake  time  constant  parameter  l~  for  EPANet 

Contaminant  traveling  as  a  slug  will  continually  approach  fresh,  uncon¬ 
taminated  pipe.  Thus,  the  rate  of  uptake  can  be  estimated  by  looking  at  the 
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initial  slope  of  the  physiochemical  uptake  model  (Equation  2).  In  other 
words,  the  derivative  of  the  physicochemical  uptake  model  is  taken  with 
respect  to  time  and  evaluating  it  at  time,  t=o,  as  shown  in  Equation  5. 


Rate  = 


Equation  5 


2.2  Reaction  rates  with  species  in  chlorinated  water 

This  problem  reduces  to  predicting  two  constants  ab  initio— the  activation 
energy  (Ea)  and  the  Arrhenius  pre-factor  (A).  The  classical  relation  be¬ 
tween  two  reactants  is  rate  =  r[A]m[B]n  where  [A]  and  [B]  are  the  molar 
concentration  of  the  two  reactants,  m  and  n  are  unknown  constant  expo- 

Ea 

nents,  and  r  =  de  kT  with  k  being  the  Boltzmann  constant  and  T  being 
temperature. 

When  a  reaction  has  rate  constants  that  obey  the  Arrhenius  equation,  both 
A  and  Ea  are  read  from  a  graph  of  ln(r)  vs.  (l/T)  at  constant  pressure,  as 
shown  in  Equation  6. 


E=-k 


f  8  In  r  ' 


5(1  IT) 


Jp 


Equation  6 


A  graph  of  this  function  immediately  yields  Ea  and  A.  This  is  used  to  vali¬ 
date  predictions  of  both  constants  using  software  such  as  Gaussian. 

Ab  initio  quantum  chemical  methods  were  used  to  investigate  the  hydroly¬ 
sis  of  an  organophosphate  (OP)  under  alkaline  conditions.  The  representa¬ 
tive  OP  consisted  of  42  atoms  arranged  in  a  complex  branched  structure. 
The  structure  of  the  OP  was  optimized  in  vacuo  to  find  the  minimum  en¬ 
ergy  configuration  using  the  6-3i+G(2d,2p)  basis  set  (Figure  2).  This  basis 
set  had  been  shown  in  previous  studies  (Wright  et  al.  2005)  to  yield  results 
that  were  in  good  agreement  with  available  experiments. 
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Figure  2.  Optimized  structure  of  an  organophosphate. 


Initial  investigations  into  possible  mechanisms  for  the  alkaline  hydrolysis 
of  OP  were  carried  out  in  the  gas  phase,  to  verily  possible  intermediate 
structures  that  had  been  obtained  from  the  scientific  literature.  Using 
these  structures  as  starting  points,  a  detailed  investigation  of  the  behavior 
of  the  OP-hydroxyl  interactions  was  undertaken  using  the  Polarizable 
Continuum  Model  (PCM)  to  simulate  liquid  water  implicitly  as  a 
continuum  with  a  known  dielectric  constant.  Transition  State  Theory 
(TST)  was  used  to  locate  saddle  points  on  the  potential  energy  surfaces 
that  correspond  to  transition  states  and  other  intermediates  in  the  alkaline 
hydrolysis  of  OP.  Pathways  were  carefully  located  and  traced  that  linked 
transition  states  to  reactants  and  products,  to  verify  that  the  mechanism 
proposed  was  physically  realizable.  Activation  energies  for  one  of  the 
possible  hydrolysis  pathways  were  obtained,  and  kinetic  equations  were 
solved  to  provide  rate  constants  for  the  reactions. 

2.3  Laboratory  validation  experiments 

2.3.1  Prophos  selection  and  analysis 

Prophos  was  selected  as  a  model  contaminant  because  it  is  hydrophobic 
and  represents  a  larger  class  of  potential  chemical  contaminants  of  con¬ 
cern.  In  addition,  it  is  similar  in  structure  to  the  organophosphates  whose 
uptake  on  surfaces  has  been  studied  by  using  molecular  modeling,  which 
allows  comparisons  between  those  results  and  the  results  of  our  laboratory 
experiments. 

Prophos  concentrations  in  solution  were  measured  at  the  University  of  Il¬ 
linois  at  Urbana-Champaign  (UIUC)  by  using  liquid  chromatog¬ 
raphy/mass  spectrometry  (LC/MS)  using  a  C-i8  column.  A  two- 
component  eluent  mixture  of  (a)  95%  water/5%  acetonitrile  and  (b)  95% 
acetonitrile/5%  water  was  run  over  a  gradient  shifting  from  25%-yo% 
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over  the  course  of  20  min.  The  LC/MS  response  for  a  set  of  standards  was 
linear  in  the  range  of  interest,  from  5  mg/L  to  750  mg/L. 

2.3.2  Reactor  vessels 

Four  types  of  pipe  reactor  vessels  were  fabricated  for  these  experiments 
from  lengths  of  piping  made  from  polyvinyl  chloride  (PVC),  copper,  and 
concrete-lined  ductile  iron  (CLDI)  as  described  below. 

Sections  of  PVC  pipe  were  cut  to  6  in.  in  length.  A  standard  PVC  primer 
was  applied  to  one  end  of  the  pipe  where  the  vessel  would  be  capped.  PVC 
cement  was  then  applied  over  the  primer.  Stock  PVC  caps  were  attached  to 
the  pipe  segments.  Care  was  taken  to  ensure  that  no  primer  or  cement  mi¬ 
grated  into  the  reactor  chamber  during  assembly. 

The  copper  reactor  vessels  were  also  prepared  using  6-in.  lengths  of  pipe. 
Stock  copper  end  caps  were  attached  using  standard  flux  and  lead-free 
solder.  Care  was  taken  to  ensure  that  no  flux  or  solder  migrated  into  the 
reactor  chamber  when  the  cap  was  attached.  Both  PVC  and  copper  reactor 
vessels  were  set  in  tip-resistant  PVC  cups  as  shown  in  Figure  3. 


Figure  3.  Reactor  vessel  (PVC)  placed  in  tip-resistant  PVC  cup. 

Two  different  types  of  CLDI  reactor  vessels  were  made.  The  first  vessel  was 
coated  with  asphaltic  sealer,  which  is  the  current  standard  for  coating  new 
pipe  materials  of  this  type.  The  second  vessel  was  unsealed  CLDI,  to  repre¬ 
sent  the  large  inventory  of  pipe  currently  in  place  in  older  water  distribu¬ 
tion  systems.  Both  vessels  were  cut  to  approximately  6-in.  length.  The  ex- 
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posed  ends  of  the  pipe  segments  were  passivated  with  marine-grade 
epoxy,  and  one  end  was  attached  to  a  standard  glass  Petri  dish.  The  labor¬ 
atory-grade  glass  in  the  Petri  dish  was  tested  and  shown  not  to  interact 
significantly  with  the  substances  of  interest.  Care  was  taken  during  the  as¬ 
sembly  to  prevent  the  epoxy  from  migrating  into  the  reactor  chamber.  The 
Petri  dish  was  of  large  enough  diameter  that  no  additional  tip-resistant 
cup  was  required  for  stability  in  the  laboratory.  Figure  4  shows  the  three 
PS  types  as  assembled. 


Figure  4.  Three  pipe  reactor  vessels  as  assembled. 


In  order  to  condition  the  pipe  materials,  deionized  water  was  placed  in  the 
vessels  and  covered  with  parafilm.  The  water  was  changed  every  1-2  days 
for  1  week.  The  PS  were  then  allowed  to  dry  overnight. 

2.3.3  Experimental  set-up 

The  prophos  uptake  experiments  were  run  at  three  different  concentra¬ 
tions:  750  mg/L,  250  mg/L,  and  125  mg/L.  Four  replicate  set-ups  of  each 
pipe  material  and  concentration  combination  were  run. 

The  prophos  solutions  were  added  to  the  pipe  segments,  after  which  the 
initial  set  of  samples  was  collected  by  removing  1  mL  of  solution  and  trans¬ 
ferring  it  to  a  vial  for  later  analysis.  The  vessels  were  then  covered  in 
parafilm  to  inhibit  contamination  and  evaporation  between  sampling 
times.  Subsequent  sampling  was  performed  periodically  over  the  next  2  to 
12  days,  with  four  replicate  set-ups  of  each  pipe  material  and  concentra¬ 
tion  combination. 
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2.4  Biofilm  studies 

Bacterial  biofilms  are  dynamic  structures  that  are  composed  primarily  of 
polysaccharides,  microbial  cells,  and  water.  Planktonic  bacteria  are  known 
to  readily  incorporate  into  biofilms,  which  can  conversely  increase  their 
resistance  to  free  chlorine  up  to  3000-fold.  The  two-fold  purpose  of  this 
task  was  to:  (1)  characterize  a  biofilm  in  relation  to  varying  free  available 
chlorine  in  a  simulated  water  distribution  system;  and  (2)  investigate  Ba¬ 
cillus  globigii  spore  uptake  within  established  biofilms. 

2.4.1  Bacillus  globigii  preparation  and  enumeration 

Dehydrated,  non-weaponized  Bacillus  globigii  spores  were  obtained  from 
Dugway  Proving  Ground  (Dugway,  UT).  Spore  viability  was  analyzed  by 
10X  serial  dilutions  over  seven  orders  of  magnitude,  followed  by  heat- 
shocking  at  80  °C  for  10  min,  and  direct  plating  100  pL  of  each  dilution 
series  onto  tryptic  soy  agar  (TSA;  Fluka,  Milwaukee,  WI).  All  dilutions 
were  plated  in  triplicate  and  incubated  for  48  hr  at  35  °C.  Orange  colony 
forming  units  (CFUs)  were  enumerated,  and  plates  that  yielded  50-200 
CFUs  were  used  in  conjunction  with  the  dilution  factors  to  calculate  the 
spore  concentration  of  the  undiluted  sample. 

2.4.2  Experimental  drinking  water  system 

Three  separate  PVC  loop  systems  were  operated  at  different  free  chlorine 
concentrations  —  0.0,  0.5,  and  1.0  mg/L  as  Cl2.  Each  pipe  loop  was  fabri¬ 
cated  out  of  schedule  40  PVC  pipe  (2.5-cm  diameter)  with  a  total  loop 
length  of  10  m.  The  total  volumetric  capacity  of  the  pipe  loop  system,  in¬ 
cluding  the  bulk  water  reservoir,  was  17  L.  Each  loop  contained  five  re¬ 
movable  20-cm  pipe  segments  for  sampling  of  biofilm  (Figure  5).  Flow  was 
controlled  by  a  magnetic  drive  centrifugal  pump  to  achieve  an  average  ve¬ 
locity  of  30  cm/sec,  measured  using  an  inline  flow  meter.  Clear  PVC  sight 
glasses  measuring  30-cm  in  length  were  installed  in  front  of  each  pump 
head  to  visually  observe  biofilm  formation  on  the  pipe  surface.  20-L  car¬ 
boy  reservoirs  were  covered  with  black  plastic  to  minimize  phototrophic 
growth.  Each  loop  system  was  operated  at  ambient  room  temperature,  ap¬ 
proximately  25  °C.  Flow  rate,  temperature,  pH,  and  free  chlorine  in  each 
pipe  loop  were  monitored  daily. 
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(A) 


Figure  5.  Simulated  PVC  water  distribution  system:  (A)  Illustration  of  pipe  loop  system 
constructed  with  standard  2.5-cm  PVC  pipe  and  fittings.  Three  separate  systems  were 
constructed,  each  having  a  total  length  of  10  m,  bulk  phase  volume  of  17  L,  operation 
temperature  of  approximately  25  °C,  and  a  flow  rate  of  30  cm/sec;  (B)  The  loop  systems  in 
operation— each  was  allowed  to  acclimate  for  63  days  as  one  unit  then  individually  poised  at 
0.0,  0.5,  and  1.0  mg/L  free  chlorine  with  sodium  hypochlorite  3  days  prior  to  the  addition  of 

B.  g/obigii  spores. 
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Prior  to  experiments,  pipe  loops  were  flushed  with  tap  water  for  a  two- 
week  period  to  remove  residual  organics  that  were  present  from  the  pipe 
manufacturing  process  and  the  joint  cement.  After  the  initial  flushing, 
each  loop  was  completely  drained  and  filled  with  dechlorinated  tap  water. 
Chlorine  was  removed  by  the  addition  of  sodium  thiosulfate  to  a  final  con¬ 
centration  of  8  mg/L.  All  three  loop  system  reservoirs  were  temporarily 
connected  to  each  other  with  a  diverter  valve,  to  allow  water  to  flow  freely 
between  each  system  for  promoting  homogeneous  biomass  formation. 
Visual  inspections  of  biofilm  formation  were  made  periodically  through 
the  sight  glasses  and  optical  densities  at  600  nm  were  taken  from  the  re¬ 
circulating  tanks  as  a  measurement  of  planktonic  growth.  Synthetic  tap 
water,  periodically  amended  with  50  mg  humic  acid,  was  used  to  stimulate 
biofilm  growth. 

After  an  additional  28  days  of  continual  operation,  the  pipe  loops  were 
drained  and  flushed  with  dechlorinated  tap  water  to  remove  residual  nu¬ 
trients  and  excess  planktonic  bacteria.  Each  pipe  loop  was  then  isolated, 
and  two  of  the  three  loops  were  amended  with  a  10%  sodium  hypochlorite 
solution  to  final  free  chlorine  concentrations  of  0.5  and  1.0  mg/L.  After 
three  additional  days  of  circulation,  each  loop  system  was  inoculated  with 
B.  globigii  to  a  final  viable  density  of  approximately  1.8  x  icH  CFUs/mL. 

2.4.3  Bulk  water  and  biofilm  sampling 

Bulk  water  and  biofilm  samples  were  taken  from  each  pipe  loop  after 
10  min  and  then  3,  7, 14,  and  21  days  of  circulation.  Liquid  samples  were 
collected  from  each  reservoir  tank  using  sterile  50-mL  pipettes.  Biofilm 
samples  were  obtained  at  predetermined  time  intervals  by  removing  one 
sampling  segment  per  loop  in  a  sequential  order.  Segments  were  replaced 
with  sterile  PVC  segments  of  the  same  dimension.  The  removed  pipe  seg¬ 
ments  were  gently  rinsed  with  sterile  ultra-pure  water  to  remove  residual 
planktonic  bacteria.  Three  2-cm  sections  were  cut  from  the  center  of  each 
pipe  segment  using  a  PVC  pipe-cutter  sterilized  with  90%  ethanol.  Bio¬ 
films  were  harvested  by  scraping.  Bulk  water  samples,  biofilm  scrapings, 
and  remaining  sample  port  segments  were  immediately  preserved  at 
-20  °C  until  further  analysis. 

2.4.4  Bacterial  enumerations 

Heterotrophic  plate  counts  (HPCs)  and  B.  globigii  viability  assays  were 
conducted  with  both  bulk  water  and  biofilm  scrapings  by  using  plating 
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methods  described  in  Section  2.3.1  with  plates  containing  plate  count  agar 
(Fluka)  and  TSA,  respectively.  Plates  were  allowed  to  incubate  for  48  hr  at 
35  °C,  and  CFUs  were  counted  on  plates  within  a  given  dilution  series  hav¬ 
ing  50-200  colonies.  Levels  of  inactivation  were  evaluated  by  plotting  the 
normalized  fraction  of  survivors  (Nt/N0)  against  exposure  time  (days). 
Values  were  expressed  as  the  total  number  of  heterotrophs  or  spores  in 
each  pipe  loop  system. 

2.4.5  Water  quality  and  biofilm  characterization 

Bulk  water-free  chlorine  concentrations  were  monitored  using  U.S.  Envi¬ 
ronmental  Protection  Agency  (USEPA)-approved  AT,AT-diethyl-p- 
phenylenediamine  (DPD)  photometric  method  (American  Public  Health 
Association  1995)  using  a  DR/2000  spectrophotometer  (HACH,  Loveland, 
CO).  Total  polysaccharide  content  of  the  biofilm  scrapings  was  determined 
with  a  modified  colorimetric  phenol-sulfuric  acid  method  (Dubois  et  al. 
1956).  Total  extra-cellular  enzyme  content  was  evaluated  using  a 
Coomassie  Plus  Better  Bradford  Protein  Assay  Kit  (Pierce  Distribution, 
Rockford,  IL)  according  to  the  manufactures  micro  test  tube  protocol.  Free 
nucleic  acids  (ssDNA,  dsDNA,  and  RNA)  present  were  quantified  using  a 
Quant-iT™  ssDNA  Assay  Kit  (Invitrogen,  Eugene,  OR)  according  to  the 
manufactures  protocol.  Bulk  water  and  biofilm  data  were  reported  as  the 
means  ±  standard  deviation  (SD)  where  n  =  3. 

2.4.6  Microscopy 

Nonviable  spores  within  biofilms  were  qualitatively  assessed  by  light  mi¬ 
croscopy  using  a  Nikon  Eclipse  E400  light  microscope  (Melville,  NY) 
equipped  with  an  Insight  digital  camera  and  Spot  imaging  software  (Diag¬ 
nostic  Instruments,  Inc.,  Sterling  Heights,  MI).  Spores  were  differentially 
stained  using  a  malachite  green  method  as  described  by  Schaeffer  and  Ful¬ 
ton  (1993)- 
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3  Results  and  Discussion 

3.1  Atomic  resolution  modeling  -  two  phases 

Our  efforts  to  describe  the  transport  of  chemicals  and  biomolecules 
through  waterworks  had  two  phases. 

3.1.1  Phase  1:  Atomic-scale  features  and  molecule  interaction  with 
surfaces 

The  first  effort  was  to  understand  how  atomic-scale  features  of  model  pipe 
surfaces  affect  the  interaction  of  the  molecules  with  the  surfaces,  and 
therefore,  determine  the  binding  and  unbinding  rates  of  the  molecules  at  a 
given  concentration.  This  phase  required  all-atom  molecular  dynamics 
simulation  using  NAMD.  Below  we  describe  how  the  development  of  all¬ 
atom  models  for  the  simulation  of  pipe  surfaces  and  contaminant  mole¬ 
cules,  as  well  as  the  results  of  methods  used  to  characterize  their  interac¬ 
tion.  This  characterization  was  performed  through  the  creation  of  atomic- 
resolution,  three-dimension  (3-D)  potential  of  mean  force  (PMF)  maps, 
which  subsume  electrical,  van  der  Waals,  and  water- mediated  interac¬ 
tions.  In  addition  to  allowing  the  effect  atomic-scale  features  on  molecule- 
surface  binding  to  be  determined,  these  maps  also  permitted  development 
of  the  multi-scale  models  used  in  the  second  phase. 

3.1.2  Phase  2:  Transport  of  molecules  near  surfaces 

In  the  second  phase,  we  investigated  explicit  transport  of  molecules  near 
surfaces  through  atomic-scale  simulation.  This  was  done  first  through  all¬ 
atom  molecular  dynamics  (MD)  simulation  using  NAMD;  however,  the 
computational  expense  of  all-atom  models  limits  its  use  to  nanoscale  sys¬ 
tems  and  short  timescales.  Thus,  we  have  developed  multi-scale  models 
using  PMF  maps  derived  from  all-atom  MD  simulations  that  are  much 
more  computationally  efficient  than  all-atom  MD  simulations  yet  lose  little 
of  the  atomic-scale  detail  required  to  describe  realistic  molecule-surface 
interaction.  These  multi-scale  models  extend  our  simulations  to  long  time- 
scales  and  large  spatial  scales,  in  addition  to  providing  the  ability  to  fully 
explore  effects  like  concentration  dependence  on  solute  binding.  We  ex¬ 
pect  the  results  and  the  methodologies  developed  from  this  work  will  be  of 
interest  not  only  for  those  modeling  solute  transport  through  public  wa¬ 
terworks,  but  also  in  the  field  of  nanofluidics,  where  the  prominence  of  so- 


ERDC  TR-12-16 


16 


lute-device  surface  interactions  can  lead  to  unexpected  behavior  and  even 
render  the  device  useless. 

3.1.3  Nanoscale  models  of  pipe  surfaces 

To  investigate  interactions  of  chemicals  and  biomolecules  with  pipe  sur¬ 
faces  at  atomic  resolution,  we  first  needed  a  description  of  the  pipe  surface 
at  the  same  level  of  detail.  To  this  end,  we  created  four  nanoscale  models 
of  pipe  surfaces  (Figure  6). 

The  first  model  represents  the  simplest  case:  An  infinitely  smooth,  fric¬ 
tionless  surface,  which  acts  as  a  hard  boundary  to  reflect  molecules  (Figure 
6A).  This  surface  did  not  consist  of  atoms;  instead  this  surface  was  mod¬ 
eled  by  applying  a  smooth  potential  to  all  the  atoms  in  the  simulation.  We 
have  implemented  this  in  NAMD,  the  molecular  dynamics  package  devel¬ 
oped  at  UIUC  through  the  grid-steered  molecular  dynamics  (G-SMD) 
method  developed  by  our  research  group. 

To  model  the  affects  of  surface  hydro phobicity  and  atomic-scale  roughness 
on  solute  adsorption,  we  created  three  all-atom  model  pipe  surfaces  com¬ 
posed  of  silica  (Si02).  To  produce  the  silica  membranes  used  in  our  re¬ 
search,  we  created  a  2.5  x  2.5  x  3.5  nm3  block  of  crystalline  silica  contain¬ 
ing  500  silicon  and  1000  oxygen  atoms  by  replicating  a  silica  unit  cell.  The 
model  membrane  was  annealed  through  a  regimen  of  simulations,  where  it 
was  first  simulated  at  high  temperatures  (7000  K)  and  gradually  cooled  to 
room  temperature.  Coordinates  of  the  silica  membrane  were  chosen  from 
three  points  in  this  regimen  to  produce  three  distinct  surfaces,  varying 
from  smooth  (Figure  6A)  and  hydrophobic  (Figure  6B),  to  a  surface  of  in¬ 
termediate  properties  (Figure  6C),  to  rough  and  hydrophilic  (Figure  6D). 


Figure  6.  Nanoscale  models  of  pipe  surfaces:  (A)  frictionless,  (B)  hydrophilic,  (C)  intermediate, 

and  (D)  hydrophobic 
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3.1.4  All-atom  model  of  DMMP 

To  study  the  transport  and  sorption/desorption  of  organophosphonates 
through  model  nanochannels,  we  chose  dimethyl-methylphosphonate 
(DMMP).  However,  at  the  time  of  our  study,  there  were  no  all-atom  pa¬ 
rameters  available  for  DMMP  that  were  compatible  with  the  Chemistry  at 
HARvard  Macromolecular  Mechanics1  (CHARMM)  force  field,  upon  which 
our  models  for  proteins  and  silica  rely.  Therefore,  we  adapted  the  united- 
atom  DMMP  parameters  developed  by  Vishnyakov  and  Neimark  (2004) 
for  use  with  the  CHARMM  force  field.  The  modified  parameters  were 
evaluated  through  numerous  MD  simulations  at  three  different  concentra¬ 
tions,  and  the  dynamics  of  the  modified  DMMP  model  were  found  to  be  in 
agreement  with  the  united-atom  model. 

3.1.5  Characterization  of  molecule-surface  interactions 

To  characterize  the  surfaces  and  determine  the  connection  between  atom¬ 
ic-scale  surface  roughness  and  solute  binding,  we  calculated  the  3-D  PMF 
between  the  DMMP  molecule  and  each  surface  by  using  umbrella  sam¬ 
pling  MD  simulations.  A  diagram  describing  this  calculation  is  shown  in 
Figure  7. 

The  PMF  at  a  given  position  represents  the  potential  energy  of  the  particle 
corresponding  to  the  average  force  exerted  on  the  molecule  at  that  posi¬ 
tion.  In  agreement  with  our  MD  flow  simulations,  we  found  that  silica  sur¬ 
face  B  is  the  most  attractive,  followed  by  surfaces  C,  D,  and  A,  with  the  lat¬ 
ter  being  the  frictionless  surface.  Silica  surface  B  shows  a  strong 
correspondence  between  the  atomic-scale  topography  and  DMMP  binding 
strength.  The  similarity  in  topography  between  silica  surface  B  and  C  leads 
us  to  conclude  that  in  addition  to  the  topography,  there  are  other  surface 
features  (e.g.,  local  surface  charge  density)  that  play  a  role  in  binding. 


1  CHARMM  is  a  versatile  and  widely  used  molecular  simulation  program  with  broad  application  to  many- 
particle  systems. 
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Figure  7.  PMF  calculations.  The  DMMP-surface  PMF  (Wms)  was  calculated  as  a  function  of  the 
position  of  the  DMMP  over  each  surface.  The  silica  surface  is  colored  by  the  value  of  the  PMF 
at  that  position  of  the  surface.  The  DMMP  is  most  attracted  to  the  regions  with  the  lowest 

PMF  (shown  in  purple). 

3.1.6  Model  nanochannels 

All-atom  molecular  dynamics  simulations  are  limited  to  length  scales  of 
tens  of  nanometers.  Furthermore,  to  increase  computational  efficiency  and 
accuracy,  we  imposed  periodic  boundary  conditions.  Thus  for  our  studies 
of  the  interactions  between  chemicals  or  biomolecules  and  pipe  surfaces, 
we  simulated  solutions  of  chemicals  or  biomolecules  flowing  in 
nanochannels.  While  there  is  a  disparity  between  the  height  of  a 
nanochannel  and  the  interior  of  a  pipe,  we  have  shown  that  the  length- 
scale  of  the  molecule-surface  interaction,  using  PMF,  is  short  enough  that 
the  majority  of  the  solution  in  the  nanochannel  is  in  the  bulk  phase  —  jus¬ 
tifying  the  use  of  our  model. 

We  created  atomic-scale  model  nanochannels  (Figure  8)  from  the  friction¬ 
less  and  silica  membranes  by  imposing  periodic  boundary  conditions  that 
defined  the  spacing  between  the  nanochannel  walls  (the  channel’s  height) 
and  made  the  channel  effectively  infinitely  long  and  wide.  Depending  on 
the  study,  the  volume  inside  the  nanochannel  was  filled  with  a  solution  of 
pure  water,  water  and  chemicals,  or  water  and  protein.  Typical  dimensions 
of  the  nanochannels  in  our  study  were  5.5  x  10  x  10  nm3. 
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Figure  8.  Diagram  of  solute  transport  through  a  silica  nanochannel.  The  average  water  profile 
calculated  from  the  all-atom  MD  simulations  is  shown  in  blue.  Solutes  are  transported  by  the 
flow  of  water  through  the  channel  and  spontaneously  bind  and  unbind  on  the  channel 
surface.  The  binding  constants  (k0/?and  k on)  calculated  from  all-atom  MD  simulations  were 

used  as  inputs  to  EPANet. 

3.1.7  NAMD  simulations  of  solute  transport 

To  induce  water  flow  through  a  nanochannel,  a  constant  pressure  differ¬ 
ence  was  created  by  applying  a  lateral  force  to  all  the  water  molecules  in 
the  system,  which  results  in  a  pressure  difference  across  the  system  (Equa¬ 
tion  7). 


A  P  =  NF/A 


Equation  7 


where: 


P  =  pressure 
N  =  atoms 
F  =  lateral  force 

A  =  area  perpendicular  to  the  applied  force 

We  simulated  the  flow  of  water  through  silica  nanochannel  B  at  five  differ¬ 
ent  pressure  gradients.  From  these  simulations,  we  determined  that  the 
mass  flux  of  water  through  the  nanochannel  will  scale  linearly,  with  the 
applied  pressure  gradient  within  the  entire  range  of  pressure  gradients 
tested  (spanning  three  orders  of  magnitude).  Furthermore,  the  simulated 
water  profiles  were  found  to  be  accurately  described  as  a  Poiseuille  flow 
through  an  infinite  channel.  However,  some  deviations  from  the  theoreti- 
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cal  predictions  were  observed  within  0.5  nm  of  surface  of  the  channel, 
where  the  breakdown  of  the  continuum  description  could  be  expected. 

For  the  remaining  nanochannels,  we  repeated  the  simulations  at  two  pres¬ 
sure  gradients  and  confirmed  the  same  behavior.  From  these  simulations, 
it  was  seen  that  the  transport  properties  of  the  silica  channels  were  very 
similar,  meaning  that  surface  properties  play  a  limited  role  in  water 
transport  through  the  nanochannel.  In  contrast  to  the  silica  surfaces,  the 
frictionless  surface  has  flat  water  velocity  profile  (a  “plug”  flow),  due  to  the 
“slip”  boundary  condition  at  the  surface. 

When  the  same  simulations  were  performed  with  the  water  replaced  by  a 
solution  of  DMMP,  we  found  similar  behavior,  but  with  an  increased  vis¬ 
cosity  of  the  fluid  that  caused  a  reduction  in  the  mass  flux.  The  addition  of 
DMMP  also  reduced  the  water  velocity  near  the  surface,  as  DMMP  ad¬ 
sorbed  on  the  channel  surfaces.  The  rough,  hydrophobic  surface  (Figure 
6D)  bound  almost  all  of  the  DMMP  in  solution;  the  intermediate  surface 
(Figure  6C)  bound  a  majority  of  the  DMMP;  and  the  smooth,  hydrophilic 
surface  (Figure  6B)  bound  less  than  half  the  DMMP  in  solution.  In  con¬ 
trast,  the  frictionless  surface  (Figure  6A)  bound  a  quarter  of  DMMP  in  so¬ 
lution.  From  these  simulations,  we  find  that  the  DMMP  molecules  are  at¬ 
tracted  to  hydrophobic  surfaces,  yet  atomic-scale  surface  roughness  also 
plays  a  large  role  in  solute  binding. 

3.1.8  Biomolecular  transport  in  nanochannels 

To  understand  the  process  of  protein  sorption  on  pipe  and  nanofluidic  de¬ 
vice  surfaces,  we  performed  the  all-atom  MD  simulations  of  protein 
transport  with  a  pressure-driven  flow  through  a  nanochannel.  With  these 
simulations,  we  characterized  (with  atomic  resolution)  the  adsorption  of  a 
model  protein  to  a  silica  surface.  Although  the  simulated  adsorption  of  the 
protein  was  found  to  be  nonspecific,  it  had  a  dramatic  effect  on  the  rate  of 
the  protein  transport. 

To  determine  the  relative  strength  of  the  protein-silica  interactions  in  dif¬ 
ferent  adsorbed  states,  we  simulated  flow-induced  desorption  of  the  pro¬ 
teins  from  the  silica  surface.  Our  analysis  of  the  protein  conformations  in 
the  adsorbed  states  did  not  reveal  any  simple  dependence  of  the  adsorp¬ 
tion  strength  on  the  size  and  composition  of  the  protein-silica  contact, 
suggesting  that  the  heterogeneity  of  the  silica  surface  may  be  an  important 
factor.  In  addition,  we  outlined  a  general  computational  method  for  char- 
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acterization  of  interactions  between  proteins  and  heterogeneous  inorganic 
surfaces  and  performed  a  proof-of-principle  demonstration  of  the  method 
using  the  frictionless  surface. 

3.1.9  Continuum  transport  model  for  solutes  in  nanochannels 

All-atom  MD  simulations  of  solute  binding  in  nanochannels  are  computa¬ 
tionally  demanding,  limiting  the  variations  in  simulation  parameters  that 
can  be  explored.  In  our  research,  we  sought  a  computational  method  by 
which  we  can  model  solute  binding  with  atomic  precision,  yet  have  less 
computational  costs  than  an  all-atom  MD  simulation,  so  we  can  fully  ex¬ 
plore  all  the  available  parameters  cheaply  and  quickly. 

To  this  end,  we  created  a  continuum  transport  model  for  solutes  in  a 
nanochannel.  In  this  model,  we  solve  the  Smoluchowski  equation,  which  is 
essentially  a  modified  diffusion  equation,  parameterized  with  the  calculat¬ 
ed  diffusion  constant  of  the  solute  and  the  calculated  PMF  of  solute  near 
each  surface.  The  Smoluchowski  model  allows  us  to  predict  both  the  flux 
of  DMMP  molecules  onto  the  channel  membrane  in  the  initial  phase  of  the 
simulations,  as  well  as  the  fraction  of  molecules  that  are  bound  to  the  sur¬ 
face  in  the  steady  state.  For  the  frictionless  surface,  the  Smoluchowski 
model  properly  predicts  both  the  transient  and  steady-state  behavior  of 
the  MD  simulations.  However,  the  model  breaks  down  for  the  silica  sur¬ 
faces,  because  the  silica  surface  PMF  varies  largely  over  the  surface  and 
the  continuum  model  has  no  particle  exclusion  so  that  crevices  in  the  silica 
surface  accumulate  too  many  particles. 

3.1.10  Brownian  dynamics  transport  model  for  solutes  in  nanochannels 

In  order  to  overcome  the  obstacles  we  faced  in  the  Smoluchowski  model, 
we  have  written  a  Brownian  dynamics  (BD)  simulator.  When  parameter¬ 
ized  with  the  calculated  diffusion  constant  of  the  solute,  the  PMF  for  the 
solute  near  the  surface,  and  the  calculated  PMF  between  two  solutes,  we 
can  perform  BD  simulations  that  recreate  the  results  of  our  all-atom  MD 
simulation.  In  effect,  we  can  simulate  the  transport  of  solutes  through 
nanofluidic  devices  with  atomic  resolution,  yet  with  a  fraction  of  the  com¬ 
putational  cost  of  an  MD  simulation.  Additionally,  we  can  now  explore  the 
dependence  of  binding  on  concentration  and  initial  conditions,  which  were 
too  computationally  demanding  for  MD.  Reference  detail  for  BD  appears 
in  Figure  9  and  Figure  10  (Carr  et  al.  2011a,  2011b). 
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Figure  9.  Comparison  of  Brownian  dynamics  (BD)  with  all-atom  molecular  dynamics,  (a)  The 
concentration  of  DMMP  bound  to  the  surface  of  each  nanochannel  system  is  a  function  of 
time,  wherein  MD  results  for  Surfaces  A,  B,  C,  and  D  (Figure  6)  are  shown  as  gray  triangles, 
red  diamonds,  black  squares,  and  blue  circles,  respectively,  with  averages  over  25  unique  BD 
trajectories  for  each  nanochannel  shown  with  lines  of  corresponding  colors;  (b)  DMMP  bound 
to  Surface  D  as  a  function  of  position— MD  results  are  on  the  left,  while  BD  results  are  on  the 
right.  The  BD  captures  the  surface  characteristics  with  atomic  resolution  detail. 


Surface  C  Surface  D 


Figure  10.  Atomic  resolution  BD  simulation  of  DMMP  advection  through  a  100  nm-long 
nanochannel  for  Surfaces  C  and  D  of  Figure  6.  A  10-mM  concentration  of  DMMP  is 
maintained  on  the  left  end  of  the  channel,  while  the  right  end  has  an  absorbing  boundary. 
DMMP  binds  more  strongly  to  the  channel  lined  with  Surface  C  than  with  Surface  D;  therefore, 
more  time  is  required  for  DMMP  to  reach  the  end  of  the  channel  in  the  first  case  than  for  the 
second  case,  despite  the  identical  water  flow  profiles. 


3.2  Laboratory  validation  experiments 

We  now  show  that  previous  research  that  attempted  to  predict  sorption  of 
contaminant  onto  pipe  material  requires  understanding  the  number  of 
contact  points  between  the  contaminant  molecule  and  the  pipe  wall.  The 
following  hypothesis  will  be  shown  to  be  inadequate.  Suppose  that  chemi¬ 
cal  agent  uptake  is  a  function  of  the  following  two  parameters: 
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•  Molecular  Weight:  This  makes  intuitive  sense  in  that  sufficiently 
large  molecules  would  have  difficulty  diffusing  into  the  molecular 
matrix  of  the  pipe  wall  material. 

•  Hydrophobicity:  This  also  makes  some  intuitive  sense  in  that  hy¬ 
drophilic  materials  were  unlikely  to  have  any  affinity  to  pipe  wall 
material.  As  a  quantitative  measure  of  hydrophobicity,  we  selected 
the  water-octanol  partition  coefficient. 

Both  of  these  measures  allowed  us  to  select  a  high-quality  simulant  (Simu¬ 
lant  P)  for  an  agent  under  study  (Agent  V).  Representative  results  from  the 
current  project  are  summarized  in  Table  l. 

Table  1.  Representative  results  for  high-quality  simulant  and  studied  agent. 

Agent  V  Simulant  P  Exp.  Error 
Molecular  Weight  267  242 

Water  octanol  partition  coefficient  123  140 

(3  w/ CLDI  w/ sealer  45%  100%  5% 

T  w /  CLDI  w/  sealer  43  hours  77  hours  2  hours 

(3  w/  cPVC  pipe  42%  4.60%  5% 

T  w/  cPVC  pipe  104  hours  107  hours  2  hours 

Here  we  see  that  the  prediction  of  matching  molecular  weight  and  hydro¬ 
phobicity  works  reasonably  well.  The  disagreement  of  measurements  be¬ 
tween  Agent  V  and  Simulant  P  are  easily  explained  by  the  sensitivity  of 
partition  coefficient  to  small  changes  in  chemical  potential.  However,  the 
yellow  line  in  Table  l  shows  that  some  difference  between  the  agent  and  its 
simulant  has  escaped  our  attention.  The  factor  of  ten  separating  the  two 
measurements  is  too  large  to  ignore. 

We  next  included  a  conjecture  that  since  chlorinated  polyvinyl  chloride 
(cPVC)  is  quite  hydrophobic  as  a  material,  there  would  be  increased  van 
der  Waals  forces  between  the  agent  molecule  and  the  pipe-wall  surface.  If 
Agent  V  is  making  contact  with  the  wall  at  two  points,  where  Simulant  P 
makes  contact  at  only  one,  we  would  expect  that  Agent  V  would  undergo 
much  more  aggressive  uptake  and  the  increase  in  chemical  potential 
would  be  somewhere  in  the  range  of  1-3  kcal/mol. 

We  then  calculated  the  magnitude  change  in  chemical  potential  and  de¬ 
termined  if  the  value  is  within  the  predicted  range.  The  experimental  re¬ 
sults  already  gave  percentage-uptake  of  contaminant  on  the  pipe  material; 
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therefore,  we  can  calculate  the  partition  coefficient  showing  the  relative 
affinity  of  the  agent  for  water  versus  pipe  material  (Equation  8). 
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^  PIPE 

r 

^  WATER 


P 

\-(3 


Equation  8 


where: 


C  =  concentration  of  the  agent  in  the  subscripted  material 
1 6  =  asymptotic  fractional  uptake  of  material  in  the  pipe 
E  =  the  chemical  potential 
k  =  Boltzmann’s  constant 
T  =  ambient  temperature  of  the  experiment 


Thus  the  change  in  chemical  potential  when  we  substitute  the  simulant  for 
the  agent  is  represented  by  Equation  g. 


AE  =  -kT\n\ 


Equation  9 


Substituting  appropriate  numbers  into  this  equation  we  get  A E  =  1.31  kcal 
per  mole.  This  is  easily  within  the  range  predicted  by  the  “one  versus  two 
contact  point”  conjecture  outlined  above.  Subsequent  chemical  modeling 
obtained  by  Ginsberg  (ERDC)  and  Hurley  (ARL)  revealed  that  this  conjec¬ 
ture  is  quite  reasonable,  as  shown  in  Figure  11  and  Figure  12  (Ginsberg  et 
al.  2008). 

In  Figure  11  and  Figure  12,  the  molecules  of  Simulant  P  and  Agent  V  are 
depicted  in  3-D  form  as  ball-and-stick  diagrams  surrounded  by  translu¬ 
cent  electron  clouds.  The  coloring  of  the  electron  cloud  shows  the  polarity 
of  the  electron  cloud  as  a  function  of  position.  Intense  red  or  intense  blue 
depict  those  areas  where  the  molecule  has  a  high  polarity  or  would  be  hy¬ 
drophilic.  The  green  areas  are  where  the  molecule  is  nonpolar  or  hydro- 
phobic. 


ERDC  TR-12-16 


25 


Figure  11.  The  3-D  depiction  of  Simulant  P.  View  at  left  shows  the  molecule  from  the  double- 
bonded  oxygen  side  (easily  seen  by  the  intense  red  spot),  and  the  right  view  is  the  obverse 
side.  Polar  regions  are  shown  in  red  and  blue;  nonpolar  regions  are  shown  in  green.  The 
molecule  is  pinwheel  shaped  and  is  likely  to  have  a  single  contact  point  with  a  nonpolar 

surface. 


Figure  12.  The  3-D  depiction  of  Agent  V.  Left  view  shows  the  molecule  from  the  double- 
bonded  oxygen  side  of  the  molecule  (easily  seen  in  the  intense  red  spot);  the  view  at  right  is 
the  obverse  side.  Polar  regions  are  shown  in  red  and  blue;  nonpolar  regions  are  shown  in 
green.  Notice  that  the  nonpolar  region  at  the  bottom  of  the  molecule  is  likely  to  have  two 

contact  points. 
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Both  molecules  have  a  double-bonded  oxygen  site  which  results  in  a  highly 
polar  red  spot.  The  molecule  is  very  unlikely  to  make  initial  contact  with 
the  pipe  wall  at  this  point.  The  Simulant  P  is  a  largely  pinwheel-shaped 
molecule  and  quite  unlikely  to  make  contact  with  the  pipe  wall  at  more 
than  one  point.  However,  Agent  V  does  have  a  large  region  opposite  the 
polar  oxygen  double  bond  that  is  quite  likely  to  make  contact  with  the  pipe 
wall  at  two  points. 

3.3  Hydrolysis  reaction  rates  and  activation  energy 

In  alkaline  solutions,  the  hydrolysis  of  OP  can  occur  via  attack  opposite  the 
ethoxide  or  thiolate  ligand.  For  the  work  described  in  this  report,  calcula¬ 
tions  were  carried  out  for  hydroxide  attack  opposite  the  thiolate  ligand.  As 
shown  in  Figure  13,  the  hydrolysis  proceeds  step-wise. 


Figure  13.  Clockwise  from  upper  left:  the  step-wise  hydrolysis  of  an  organophosphonate  via 
attack  opposite  the  thiolate  ligand.  (The  representative  colors  are  yellow  =  phosphorus,  dark 
yellow  =  sulfur;  red  =  oxygen,  blue  =  nitrogen,  grey  =  carbon,  and  white  =  hydrogen.) 


A  hydroxide  molecule  must  first  align  itself  with  the  phosphonate  center  in 
the  OP  molecule.  It  subsequently  binds  to  the  phosphorus  atom,  destabi¬ 
lizing  the  phosphorus-sulfur  (P-S)  bond  and  resulting  in  bond  cleavage. 
The  rate-determining-step  (i.e.,  the  slowest  of  the  reactions)  for  the  pro¬ 
posed  mechanisms  is  the  initial  alignment  of  the  hydroxyl  radical  with  the 
phosphonate  center  of  the  OP  molecule.  The  transition  of  OP  to  a  bound 
state  with  hydroxide  through  the  aligned  intermediate  A  requires  activa- 
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tion  energy  of  21.08  kcal/mol.  Subsequent  steps  in  the  hydrolysis  pathway 
are  energetically  favored  (Figure  14). 


Figure  14.  Changes  in  the  energy  of  intermediate  states  in  the  hydrolysis  pathway. 

Reaction  rate  constants  were  also  calculated  for  each  step  of  the  hydrolysis 
reaction  shown  in  Figure  14.  The  values  were  calculated  by  using  TST  ac¬ 
cording  to  Equation  10. 


kT 

rv(T)  =  (—)e 
h 


-GTs  Uj) 
RT 


Equation  10 


where: 


rij(T)  =  rate  constant  of  the  transition  i  to  j 
k  =  Boltzmann  constant 
h  =  Planck’s  constant 
R  =  gas  constant 

T  =  temperature,  298.15  K  for  these  calculations 

The  calculated  rate  constants  are  summarized  in  Table  2.  As  expected,  the 
rate  constant  for  the  alignment  of  a  hydroxide  molecule  with  the  OP  is  sig¬ 
nificantly  slower  (by  several  orders  of  magnitude)  than  the  other  reactions. 
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Table  2.  Reaction  rate  constants  for  OP  hydrolysis  via  hydroxide  attack 
opposite  the  thiolate  ligand. 


Reaction 

Rate  Constant  [sec1] 

V  ->  A 

2.2  x  1013 

A  ->  B 

7.7  x  1035 

B  ->  C 

3.2  x  1014 

C  ->  D 

4.3  x  1028 

Calculations  of  hydrolysis  activation  energy  and  reaction  rate  constants 
under  varying  pH  and  temperatures  were  performed  for  comparison, 
based  on  previously  reported  laboratory  data  (Morrissey  et  al.  2009).  Half- 
lives  (tv2)  at  pH  8  and  14  were  obtained  from  the  report  and  used  to  calcu¬ 
late  rate  constants  according  to  Equation  11,  assuming  a  first-order  reac¬ 
tion. 


In  2 
r  = - 


Equation  11 


Reaction  rate  constants  were  also  calculated  from  values  for  the  time  re¬ 
quired  to  achieve  99.9%  hydrolysis  at  different  temperatures  according  to 
Equation  12. 


C  =  Coert 


Equation  12 


where: 


C  =  the  concentration  at  time  t 
Co  =  the  initial  concentration 
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Table  3  lists  the  first-order  reaction  rate  constants  obtained  for  OP  hydrol¬ 
ysis  from  the  above  calculations. 


Table  3.  Reaction  rate  constants  for  OP  hydrolysis  via  hydroxide  attack 
opposite  the  thiolate  ligand. 


pH 

Temperature  [K] 

Rate  Constant  [sec1] 

8* 

298 

1.1  x  10-6 

14* 

298 

9.6  x  10-3 

8A 

278 

6.6  x  10'8 

8A 

288 

2.7  x  10-7 

8A 

298 

1.1  X  IQ'6 

*  Data  taken  from  Morrissey  et  al.  1999a; A  Data  taken  from  Morrissey  et  al.  1999b 

These  reaction  rate  constants  are  considerably  more  rapid  than  those 
obtained  from  the  computational  predictions  for  the  rate  limiting  step. 
However,  two  key  points  must  be  made.  Firstly,  the  calculations  for  the 
modeling  focused  only  on  one  hydroxide  attack  site.  Hydroxide  attack  may 
also  align  and  subsequently  bind  to  the  phosphorus  opposite  the  ethoxide 
ligand.  Additionally,  the  modeling  calculations  were  in  essentially 
infinitely  basic  solution.  As  demonstrated  by  the  laboratory  data,  the  pH 
(and  thus  the  availability  of  hydroxide)  plays  a  significant  role  in  reaction 
rate  constants. 

In  summary,  a  direct  comparison  of  reaction  rate  constants  is  difficult  due 
to  the  variations  in  conditions  assumed  for  the  calculations.  However, 
activation  energy  for  hydrolysis  is  less  dependent  on  the  conditions  of  the 
system.  For  the  laboratory  data,  activation  energy  can  be  calculated  from 
reaction  rate  constants  determined  at  different  temperatures  as  shown  by 
Equation  13  and  Equation  14. 


r—Ae  RT 


Equation  13 


In  r 


E 

— —  +  \nA 

RT 


Equation  14 


where: 


r  =  amount  of  reagent  present 
Ea  =  activation  energy 
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R  =  Boltzman’s  constant 
T  =  time 
A  =  Constant 


So,  a  plot  of  In  r  vs.  l/T  gives  a  slope  of  -  . 

Multiplying  by  -R  yields  the  activation  energy.  For  the  hydrolysis  of  an  or- 
ganophosphate  data  we  have  available,  the  activation  energy  is 
23.06  kcal/mol,  which  is  within  10%  of  the  value  obtained  from  the  com¬ 
puter  simulations. 

3.4  Biofilm  characterization  and  spore  fate 

The  biofilm  environments  were  characterized  by  daily  measuring  of  free 
chlorine  concentration,  temperature,  and  pH  during  the  21-day  period. 
Measured  values  from  the  days  on  which  bacterial  enumerations  took 
place  are  reported  in  Table  4.  Free  chlorine  values  remained  relatively  sta¬ 
ble  throughout  the  study.  Concentrations  in  the  0.5  mg/L  and  the  1.0 
mg/L  amended  loops  varied  throughout  the  study  by  0.1/L  and  0.3  mg/L, 
respectively.  No  free  chlorine  was  detected  in  the  untreated  pipe  loop  until 
after  the  addition  of  sodium  hypochlorite  on  Day  14.  Temperatures  fluctu¬ 
ated  less  than  0.5  °C  between  each  of  the  three  pipe  loops  on  a  given  day. 
Slight  variations  in  pH  were  observed  among  the  pipe  loops;  however, 
readings  remain  well  within  the  range  of  a  typical  potable  water  distribu¬ 
tion  system  (World  Health  Organization  2007).  Due  to  the  tremendous 
physiochemical  tolerances  associated  with  Bacillus  spores  (Nicholson  et  al. 
2000),  the  slight  variability  in  these  physiochemical  parameters  were  not 
believed  to  have  a  significant  effect  on  B.  globigii  resistance  within  any 
given  pipe  loop  over  the  course  of  the  study.  Turbidity  was  monitored  pri¬ 
or  to  the  start  of  the  study  as  an  indicator  of  biological  activity.  Each  of  the 
pipe  loop  systems  increased  in  absorbance  from  0.01  to  0.03  after  three 
weeks  of  circulation  in  the  presence  of  humic  acid.  Consistency  of  these 
ratios  indicated  that  each  pipe  loop  contained  similar  levels  of  biological 
activity.  This  consistency  was  confirmed  by  enumerating  the  bulk  water 
heterotrophs,  which  revealed  that  each  system  contained  approximately 
3.4  x  io5CFU/mL  prior  to  flushing  the  system.  The  untreated  pipe  loop 
harbored  1.4  x  icH  CFU/cm2  at  the  start  of  the  study,  which  was  indicative 
of  a  mature  water  distribution  system  (Schwartz  et  al.  2003).  These  data 
demonstrated  that  the  conditions  in  each  of  the  pipe  loops  were  similar  at 
the  onset  of  the  study  and  that  each  remained  stable  throughout  the  21- 
day  circulation  period. 
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Table  4.  PVC  pipe  loops,  bulk  water  physiological  conditions. 


Time  (days)b 

0  3  7  14  21 


0.0  mg/L  free  chlorine'1 


Free  Cl'  (mg/L) 

0.0 

0.0 

0.0 

0.0 

0.9 

Temperature  (°C) 

23.3 

25.4 

22.8 

26.4 

23.2 

pH 

7.5 

7.2 

7.4 

7.4 

7.6 

0.5  mg/L  free  chlorine'1 

Free  CF  (mg/L) 

0.5 

0.4 

0.4 

0.5 

1.0 

Temperature  (°C) 

23.6 

25.9 

23.1 

26.9 

23.4 

pH 

7.7 

7.4 

7.6 

7.6 

7.7 

1.0  mg/L  free  chlorine 

Free  CF  (mg/L) 

1.0 

1.1 

0.9 

1.1 

1.2 

Temperature  (°C) 

23.2 

25.5 

23.0 

26.8 

23.3 

pH 

7.7 

7.8 

7.9 

7.9 

7.9 

“Pipe  loop  systems  were  amended  to  a  target  concentration  of 
1.0  mg/L  free  chlorine  on  day  14. 

bDays  after  the  addition  of  3.0  x  10s  B.  globigii  spores  to  each 
PVC  pipe  loop  system  (Time  0=10  min). 

The  biofilms  were  characterized  by  measuring  total  polysaccharide,  free 
protein,  and  extracellular  nucleic  acid  concentrations  (Table  5).  In  the 
chlorine-free,  untreated  pipe  loop,  the  polysaccharide  concentrations  in¬ 
creased  steadily  over  the  course  of  the  study  to  a  maximum  concentration 
of  3.6  pg/cm2  after  14  days  of  circulation.  On  Day  14,  free  protein  and  nu¬ 
cleic  acids  also  accumulated  to  maximum  concentrations  of  0.3  pg/cm2 
and  2.4  pg/cm2 ,  respectively.  Poising  the  system  at  1.0  mg/L  free  chlorine 
on  Day  14  decreased  the  polysaccharide  concentration  to  1.6  pg/cm2;  both 
extracellular  proteins  and  free  nucleic  acids  were  completely  oxidized  after 
seven  additional  days  of  circulation.  The  biofilms  within  the  chlorine- 
treated  pipe  loops  remained  intact,  and  no  statistical  variation  in  polysac¬ 
charide  concentrations  were  observed  over  the  course  of  the  study.  This 
finding  was  consistent  with  studies  that  have  demonstrated  that  free  chlo¬ 
rine  concentrations  less  than  5  mg/L  have  little  effect  on  biofilm  stability 
(Chen  and  Stewart  2000;  Daly  et  al.  1989).  The  nearly  two-fold  decrease  in 
the  biofilm  composition  observed  in  the  untreated  pipe  loop  after  the  addi¬ 
tion  of  sodium  hypochlorite  on  Day  14  was  likely  due  to  detachment 
caused  by  sudden  environmental  stress  (Garny  et  al.  2009). 
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Table  5.  PVC  pipe  loops,  biofilm  conditions. 


0.0  mg/L  free 
chlorinea,b 

Polysaccharide  (pg/cm2) 
Protein  (pg/cm2) 

Nucleic  acids  (ng/cm2) 

0 

3 

Time  (days)' 
7 

14 

21 

0.6  ±0.1 

0.1  ±0.0 

0.8  ±0.4 

1.2  ±0.4 
0.1  ±0.0 
1.7  ±  0.7 

1.4  ±0.5 

0.2  ±0.1 

1.1  ±0.4 

3.6  ±0.5 
0.3  ±  0.0 

2.4  ±  1.2 

1.6  ±  0.1 

0.0  ±0.0 

0.0  ±0.0 

0.5  mg/L  free 
chlorineab 

Polysaccharide  (pg/cm2) 
Protein  (pg/cm2) 

Nucleic  acids  (ng/cm2) 

0.6  ±  0.1 

0.1  ±0.0 

0.0  ±0.0 

0.6  ±0.2 

0.1  ±0.0 

0.0  ±0.0 

0.8  ±0.2 

0.1  ±0.0 

0.0  ±0.0 

0.4  ±0.2 

0.1  ±0.0 

0.0  ±0.0 

0.4  ±0.2 

0.0  ±0.0 

0.0  ±0.0 

1.0  mg/L  free  chlorine'1 

Polysaccharide  (pg/cm2) 
Protein  (pg/cm  ) 

Nucleic  acids  (ng/cm2) 

0.4  ±0.1 

0.0  ±0.0 

0.0  ±0.0 

0.4  ±0.2 

0.0  ±0.0 

0.0  ±0.0 

0.6  ±0.3 

0.0  ±0.0 

0.0  ±0.0 

0.8  ±0.2 

0.0  ±0.0 

0.0  ±0.0 

0.6  ±0.2 

0.0  ±0.0 

0.0  ±0.0 

3Data  are  means  ±  SD  (n  =  3). 

bPipe  loop  systems  were  amended  to  a  target  concentration  of  1 .0  mg/L  free  chlorine 
on  day  14. 

cDays  after  the  addition  of  4.0  x  1 0s  B.  globigii  spores  to  each  PVC  pipe  loop  system 
(Time  0=10  min). 


Free  protein  remained  steady  at  o.i  pg/cm2  in  the  0.5  mg/L  free  chlorine- 
amended  pipe  loop,  until  the  chlorine  concentration  was  increased  to  1.0 
mg/L  on  Day  14,  when  proteins  were  completely  oxidized.  No  proteins  or 
nucleic  acids  were  detected  in  the  1.0  mg/L  free  chlorine  amended  pipe 
loop  throughout  the  study.  Collectively,  the  biofilm  physiological  condi¬ 
tions  clearly  demonstrated  the  ability  of  the  PVC  pipe  biofilm  to  grow  un¬ 
der  oligotrophic  conditions.  Furthermore,  free  chlorine  concentrations  of 
0.5  mg/L  were  sufficient  to  penetrate  and  suppress  polysaccharide  for¬ 
mation  in  addition  to  oxidize  nucleic  acids;  however,  free  chlorine  concen¬ 
trations  greater  than  0.5  mg/L  were  required  to  completely  oxidize  biofilm 
proteins. 

The  concentrations  of  heterotrophic  bacteria  in  the  bulk  water  were  also 
monitored.  HPCs  were  conducted  on  the  pipe  loop  bulk  water  prior  to  the 
addition  of  sodium  hypochlorite,  and  approximately  3.4  x  ic>3  CFUs/mL 
were  observed  in  each  system  after  flushing  with  dechlorinated  tap  water. 
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Heterotrophs  remained  stable  in  the  chlorine-free,  untreated  system  until 
the  addition  of  1.0  mg/L  free  chlorine  on  Day  14  (Figure  15),  after  which 
the  total  number  of  heterotrophs  was  below  detection  after  7  additional 
days.  In  the  0.5  mg/L  amended  pipe  loop,  heterotrophs  decreased  1.6-logio 
within  the  3-day  stabilization  period  prior  to  inoculation.  No  viable  HPCs 
were  detected  in  the  1.0  mg/L  free  chlorine  treated  pipe  loop  throughout 
the  21-day  study,  indicating  complete  inactivation. 


Figure  15.  Susceptibility  of  heterotrophic  bacteria  and  B.  globigii spores  to  free  chlorine  in 
PVC  pipe  loop  bulk  water,  presented  as  the  logio  of  the  total  surviving  fraction.  Approximately 
3.0  x  108  viable  B.  globigii spores  were  added  to  each  pipe  loop  system  on  Day  Zero.  Free 
chlorine  in  all  systems  was  poised  at  1.0  mg/L  on  Day  14.  (Symbols:  ♦  =  heterotrophs, 
chlorine  untreated;  0  =  B.  globigii,  chlorine  untreated;  ■  =  heterotrophs,  0.5  mg/L  free 
chlorine;  □  =  B.  globigii,  0.5  mg/L  free  chlorine;  •  =  heterotrophs,  1.0  mg/L  free  chlorine;  and 
o  =  B.  globigii,  1.0  mg/L  free  chlorine.  Data  are  means  +  SD  (n=3);  some  error  bars  are  not 

visableduetosize. 


The  concentration  of  viable  B.  globigii  spores  in  the  bulk  solution  of  each 
pipe  was  also  monitored.  In  the  chlorine-free,  untreated  pipe,  the  spore 
concentration  was  1.5  x  104  ±  1.2  x  103  CFUs/mL  and  remained  steady  until 
the  introduction  of  1.0  mg/L  free  chlorine  on  Day  14  (Figure  16).  Within 
seven  days  of  spiking  the  chlorine,  spore  viability  was  reduced  to  levels  be¬ 
low  the  detection  limit.  The  spore  concentration  in  the  0.5  mg/L  pipe  loop 
was  2.3  x  104  ±  1.2  x  103  CFUs/mL  after  10  min  of  circulation.  Despite  ex- 
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posure  to  0.5  mg/L  of  chlorine  for  10  min,  no  difference  in  spore  concen¬ 
tration  was  observed  between  the  untreated  and  0.5  mg/L  free  chlorine 
amended  pipe  loops,  indicating  no  inactivation  had  occurred  in  the  initial 
timeframe.  Subsequently,  inactivation  of  greater  than  four  logs  was  ob¬ 
served  within  14  days.  In  the  pipe  loop  containing  1.0  mg/L  free  chlorine, 
no  viable  spores  were  detected  within  the  bulk  water  throughout  the  study. 


Figure  16.  Susceptibility  of  heterotrophic  bacteria  and  B.  globigii spores  to  free  chlorine  in 
PVC  pipe  biofilm,  presented  as  the  logio  of  the  total  surviving  fraction.  Approximately  3.0  x 
108  B.  globigii s po res  were  added  to  each  pipe  loop  system  on  day  zero.  Free  chlorine  in  all 
systems  was  poised  at  1.0  mg/L  on  day  14.  Symbols:  ♦  heterotrophs,  chlorine  untreated;  0  B. 
globigii,  chlorine  untreated;  ■  heterotrophs,  0.5  mg/L  free  chlorine;  □  B.  globigii,  0.5  mg/L 
free  chlorine;  •  heterotrophs,  1.0  mg/L  free  chlorine;  and  o  B.  globigii,  1.0  mg/L  free 
chlorine.  Data  are  means  ±  SD  (n=3),  some  error  bars  are  not  visable  due  to  size. 


Biofilm  HPCs  in  the  chlorine-unamended  pipe  loop  were  1.4  x  103  ±  2.2  x 
102  CFUs/cm2  at  the  start  of  the  study  and  remained  stable  until  chlorine 
treatment  (Figure  16).  Addition  of  1.0  mg/L  free  chlorine  into  the  non- 
chlorinated  pipe  loop  resulted  in  a  reduction  in  HPCs  to  non-detectable 
levels  within  7  days.  A  1.0-logio  reduction  in  heterotrophs  was  observed 
after  3  days  circulation  in  the  0.5  mg/L  free  chlorine  treated  pipe  loop,  and 
complete  inactivation  was  achieved  after  7  days  of  incubation.  Consistent 
with  published  studies,  sessile  heterotrophs  were  more  resistant  (nearly  3- 


ERDC  TR-12-16 


35 


times)  to  chlorine  inactivation  than  planktonic  forms  (Murga  et  al.  2001; 
Park  et  al.  2001;  Steed  and  Falkinham  2006;  Szabo  et  al.  2006;  DeQueiroz 
and  Day  2007;  Bressler  et  al.  2009).  As  expected,  elevated  free  chlorine 
concentrations  increased  the  rate  of  heterotroph  disinfection.  As  seen  with 
the  bulk  water,  no  CFUs  were  detected  within  the  biofilm  of  the  pipe  loop 
amended  with  1.0  mg/L  free  chlorine  throughout  the  21-day  study,  indi¬ 
cating  100%  inactivation  within  the  3  days  prior  to  spore  inoculation. 

B.  globigii  spores  were  found  to  associate  rapidly  with  the  pipe  surface 
biofilm.  After  10  min  of  circulation,  285  ±  92  and  480  ±  6  viable 
spores/cm2  of  pipe  and  0.5  mg/L  free  chlorine  loops,  respectively,  were 
detected  in  biofilm  samples  from  the  untreated  pipe  loop  (Figure  16).  After 
3  days  of  circulation,  the  number  of  spores  associated  with  the  biofilm  in 
the  untreated  pipe  loop  increased  to  1.3  x  lcri  ±  205  spores/cm2  (0.7-logio 
increase),  which  represented  roughly  3%  of  the  total  spores  added  to  the 
system.  An  aqueous  infective  dose  is  estimated  to  be  171  spores/L,  at  a 
consumption  rate  of  5  L/day  for  seven  days  (Burrows  and  Renner  1999). 
Studies  have  shown  that  more  than  90%  of  a  pipe  biofilm  can  be  sloughed 
from  the  surface  (Garny  et  al.  2009),  suggesting  that  spore  concentrations 
observed  within  the  PVC  pipe  biofilm  during  the  current  study  were  theo¬ 
retically  sufficient  to  cause  infection  in  the  event  of  biofilm  detachment. 

Spore  numbers  remained  statistically  consistent  in  the  untreated  pipe  loop 
until  the  addition  of  1.0  mg/L  free  chlorine  to  the  system  on  Day  14.  After 
chlorine  treatment,  B.  globigii  spores  were  not  detected  after  7  additional 
days  of  circulation,  which  represented  a  7.1-logio  reduction  in  the  total  ses¬ 
sile  spore  population  (Figure  16).  In  the  0.5  mg/L  free  chlorine  loop  sys¬ 
tem,  only  a  2.7-logio  reduction  in  sessile  spores  was  observed  after  14  days 
of  circulation  when  compared  to  a  decrease  below  detection  in  the  bulk 
water  phase.  This  observation  indicated  a  decrease  in  disinfection  efficien¬ 
cy  within  the  biofilm,  likely  due  to  protection  of  enmeshed  spores  from 
chlorine  (De  Beer  et  al.  1994;  Szabo  et  al.  2006). 

Interestingly,  the  biofilms  that  developed  in  both  the  0.5  mg/L  and  the  1.0 
mg/L  free  chlorine  pipe  loops  had  approximately  equal  concentrations  of 
polysaccharides,  yet  no  viable  spores  were  detected  within  the  biofilm  of 
1.0  mg/L  amended  pipe  loop  throughout  the  study.  Light  microscopy  re¬ 
vealed  spores  were  embedded  within  the  biofilm.  The  lack  of  viable  spores 
within  the  bulk  water  indicated  spores  were  likely  inactivated  prior  to  in¬ 
corporation  into  the  biofilm.  In  the  0.5  mg/L  free  chlorine  pipe  loop,  the 
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slower  inactivation  rate  allowed  viable  spores  to  integrate  within  the  bio¬ 
film,  which  resulted  in  elevated  persistence  within  the  loop  system. 

Although  1.0  mg/L  free  chlorine  was  demonstrated  to  completely  inacti¬ 
vate  1.8  x  104 B.  globigii  spores/mL  within  the  PVC  pipe  systems,  it  is  not¬ 
ed  that  these  experiments  were  conducted  at  room  temperature,  and  chlo¬ 
rine  inactivation  at  typical  water  distribution  temperatures  is  expected  to 
reduce  logio inactivation  ratios  appreciably  (Ndiongue  et  al.  2005).  Rice  et 
al.  (2005)  observed  approximately  a  three-fold  increase  in  2-  logio  inacti¬ 
vation  values  for  B.  anthracis  Sterne  and  B.  cereus  when  temperature  was 
decreased  from  23°C  to  5°C  in  the  presence  of  free  chlorine.  When  treated 
with  monochloramine,  Rose  et  al.  (2007)  reported  increased  2-logio  inacti¬ 
vation  times  of  eight-fold  for  B.  anthracis  Ames  and  seven-fold  for  B. 
anthracis  Sterne,  when  temperatures  decreased  from  25  °C  to  5  °C.  These 
reports  suggest  that  spore  inactivation  intervals  for  a  given  logio  reduction 
within  the  pipe  loop  systems  described  in  this  study  would  increase  signif¬ 
icantly  at  lower  water  temperatures. 
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4  Implementation  into  EPANet  MSX 

EPANet  is  quite  uneven  in  its  ability  to  accept  the  results  of  our  research. 
Some  capabilities  to  model  applicable  ordinary  differential  equations  are 
built  in.  For  the  reaction  rates  of  bacillus  species  in  chlorinated  water,  the 
software  merely  requires  a  supplemental  library  from  which  it  can  draw 
kinetic  parameters  based  on  the  contaminant  selected. 

For  the  sorption  /  desorption  model  developed  through  the  use  of  NAMD, 
a  general  description  of  how  to  derive  the  ordinary  differential  equation  of 
interest  is  based  on  the  form  of  the  function/is  beyond  the  scope  of  this 
document.  However,  it  is  instructive  to  quickly  summarize  how  this  is 
done  for  our  previously  measured  uptake  equation,  reproduced  here  as 
Equation  15. 


Pit) 


Equation  15 


Here  we  can  observe  that  the  system  has  two  components  (water  and  pipe) 
and  that  the  dynamics  are  first  order  (one  time  constant,  T).  Hence,  we  can 
readily  derive  the  first  order,  ordinary  differential  equation  governing  stat¬ 
ic  sorption  and  desorption  in  static  conditions,  shown  by  Equation  16. 

C{t) 

Pit) 


-p/r  (i-/?)/r 
p/r  ~(i-p)/r 


ao 

Pit) 


f  rn  1  1  o  t  i  rv  n  d  (Pt. 


Where 

C(f)  =  contaminant  concentration  in  water 
Pit)  =  the  concentration  in  the  pipe  wall 
P  =  asymptotic  uptake  and 
T  =  time  constant 

The  differential  equation  has  eigenvalues  at  o  and  -1  /T.  The  associated  ei¬ 
genvectors  are  [1-  (3,  [1]  and  [1,-1]  respectively.  This  makes  intuitive  sense 
because  the  first  eigenvalue  of  o  gives  an  eigenvector  pointing  along  the 
the  stabile  equilibrium  line  in  the  phase  space.  This  is  consistent  with 
thermodynamic  partition  and  predicted  by  the  uptake  equation.  The  se- 
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cond  eigenvalue  shows  the  correct  time  constant  of  the  system  as  it  moves 
toward  equilibrium.  The  direction  of  the  second  eigenvector  forces  the  sys¬ 
tem  to  move  along  lines  in  the  phase  space  where  C(f)+P(f)  =  constant  so 
that  the  mass  of  total  contaminant  is  conserved. 

If  we  were  writing  the  pipe  segment  hydraulic  simulation  from  scratch,  we 
would  then  proceed  to  modify  this  differential  equation  to  capture  water 
flow  and  aggregation  at  the  pipe  segment  level.  The  resulting  differential 
equation  would  then  be  expressed  as  a  difference  equation.  Having  re¬ 
viewed  the  source  code  of  EPANet  2.0,  these  changes  appear  in  the  file 
quality,  c,  in  the  subroutines  ratecoeffs,  piper  ate,  and  ivallrate. 
EPANet  uses  a  Lagrangian  Time-Drive  Method  to  calculate  system  up¬ 
dates.  Reaction  equations  and  flow  equations  are  run  in  alternating  se¬ 
quence,  with  the  spatial  size  of  pipe  lengths  determined  adaptively  by  cre¬ 
ating  new  lengths  when  the  concentration  of  a  water  constituent  in 
adjacent  pipe  lengths  exceeds  a  fixed  threshold. 

4.1  Simulations  of  fate  and  transport  using  EPANet  MSX 

The  fate  and  transport  equation  above  was  implemented  in  EPANet  MSX 
(multi-species  extension),  which  is  under  development  at  USEPA’s  Na¬ 
tional  Homeland  Research  Center  and  at  the  University  of  Cincinnati.  The 
preamble  code  is  a  direct  quote  of  the  equation  and  is  given  in  the  follow¬ 
ing  diagram.  Rate  constants  were  taken  from  ERDC-CERL  molecular 
modeling  measurements  confirmed  by  laboratory  analysis. 

As  a  result,  EPANet  MSX  was  able  to  render  animations  of  simulation  re¬ 
sults,  as  shown  in  Table  6,  to  clearly  show  the  difference  between  contam¬ 
inants  that  are  hydrophilic  (easily  soluble  in  water)  and  those  that  are  hy¬ 
drophobic  (likely  to  stick  to  the  pipe  surface). 

When  the  hydrophobic  contaminant  sticks  to  the  pipe,  the  time  constants 
describing  uptake  and  release  are  significantly  longer  than  time  constants 
describing  hydraulic  transport.  Therefore,  the  contamination  event  has  an 
expected  duration  of  weeks,  not  hours.  This  duration  difference  is  also 
summarized  in  Figure  17.. 
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Table  6.  ERDC-developed  EPANET  MSX  physicochemical  model  simulations,  showing  the  end 
state  after  24  hours  of  a  contamination  event  in  two  water  distribution  system  layouts 
comparing  hydrophilic  contaminant  distribution  in  water  and  on  pipe  walls. 


At  24  Hours 
After 

Contamination 


Hydrophilic  Contaminant 
(concentration  in  bulk  water) 


Hydrophobic  Contaminant 
(concentration  in  bulk  water) 


Hydrophobic  Contaminant 
(concentration  on  pipe) 


Subdivision 
(no  tanks) 


V, 


Mains 
(with  tanks) 


l 


Problem:  Hydrophobic  Contaminants  Undergo 
Adsorption  and  Desorption  from  Pipe  Material 


■  No  chlorine  interaction 

■  CHPPM  drinking  water  limit  at  1.2  c  Q01 
x  10'8  (  =  12  ppb) 

■  Time  to  reach  drinking  water  limit:  § 

►  No  uptake  ~  0.5  hr  °  t3 

►  With  uptake  ~  240  hr 


Concentration  in  w  ater 
(as  seen  at  target) 

- with  Uptake 

- no  Uptake 

Agent:  V  (Concentrated) 

Pipe:  Cement-Lined  Ductile  Iron 
Pipe  Inner  Diameter:  3.125" 
Pipe  Length:  3600  ft 
Flow:  1  foot  per  second 
Dose:  60  gallons 
Injection  Rate:  12  gpm 


I 


0  9)  100  150  200  250  300 

Time  [minutes] 


BUILDING  STRONG, 


Figure  17.  Hydrophilic  water  contamination  would  take  a  matter  of  hours,  not  days. 
Hydrophobic  contamination  would  linger  for  a  number  of  days. 
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5  Future  Research  Efforts 

5.1  Water  Integrated  Security  Test  Laboratory 

Future  goals  include  expansion  and  validation  of  the  results  obtained  from 
this  work,  with  increased  focus  on  developing  an  integrated  system  of  pro¬ 
tection  (Figure  18).  There  is  a  need  to  develop  reliable  systems  with  the 
ability  to  detect  and  respond  to  water  contamination  in  a  timely  manner  to 
keep  water  distribution  systems  protected  from  intentional  and  accidental 
poisoning. 

In  order  to  examine  the  reliability  of  system  components,  research  on  de¬ 
tection,  prevention,  and  mitigation  of  contamination  of  water  distribution 
systems  in  a  small  scale  will  be  performed  at  the  Water  Integrated  Security 
Test  Laboratory  (WISTL)  located  at  ERDC-CERL  Champaign,  IL.  The 
WISL  allows  replication  of  different  water  conditions  around  the  world 
and  simulation  of  possible  water  contamination  events.  Such  events  will  be 
detected  by  a  real-time  early  warning  system  (HACH  Homeland  Security 
Technologies  Guardian  Blue®)  and  a  real-time  treatment  station  with  a 
high-quality,  unique  filtration  system  based  on  Teslin®,  a  synthetic  print¬ 
ing  medium  manufactured  by  PPG  Industries. 


& 

Threat 

Agent 


Water  Facility 


Office  or  Residence 
with  Drinking  Water  Hookup 


Backflow  Attack 
Staging  Ground 


Control  Signed  "^i 


Real-Time  Decentralized 
Treatment  Doser 


Icon  Facility, 
Political  Target,  or 
Military  Target 


Figure  18.  Integrated  water  system  protection. 


ERDC  TR-12-16 


41 


5.2  Secondary  aerosol  modeling 

When  a  biological  warfare  attack  is  launched  through  aerosolization,  the 
initial  contaminant  plume  is  referred  to  as  the  primary  aerosol.  The  pri¬ 
mary  aerosol  typically  accounts  for  the  greatest  risk  of  illnesses  and  death 
associated  with  a  biological  attack.  Under  normal  circumstances,  a  biologi¬ 
cal  aerosol  can  be  dispersed  and  begin  to  be  deposited  on  surfaces  within 
hours  of  release,  depending  on  meteorological  conditions  and  spore  aero¬ 
dynamics  (Meselson  et  al.  1994).  Evidence  suggests  that  biological  agents 
may  be  resuspended  from  the  ground  or  other  surfaces  (e.g.,  HVAC  ducts) 
and  continue  to  exist  in  the  atmosphere  as  a  secondary  aerosol  (Inglesby 
et  al.  1999;  Weis  et  al.  2002).  Thus,  there  is  potential  that  persistent  sec¬ 
ondary  aerosols  will  continue  to  pose  a  significant  health  hazard  to  build¬ 
ing  occupants. 

Experimental  work  has  been  done  to  characterize  the  resuspension  of  par¬ 
ticles  from  various  surfaces  common  to  an  office  building  under  various 
flow  characteristics  (Mukai  et  al.  2009).  While  it  is  useful  to  do  experi¬ 
mental  work  such  as  this,  it  is  impossible  to  fully  understand  the  scope  of 
the  problem  due  to  many  variables  in  airflow  characteristics  as  well  as  dif¬ 
fering  rates  of  re-aerosolization  among  particles  of  interest.  Initial  theoret¬ 
ical  work  has  been  completed  to  determine  the  resuspension  term  from 
indoor  surfaces  as  applied  to  commercially  available  computational  fluid 
dynamics  (CFD)  codes  (Kim  et  al.  2010). 

While  there  is  significant  importance  in  modeling  sorption  and  desorption 
of  contaminants  from  pipe  walls  in  an  aqueous  system,  there  is  equal  im¬ 
portance  in  modeling  the  same  trends  in  air.  Chemical  and  biological  de¬ 
fense  is  a  key  focus  within  the  Army,  as  listed  as  an  Army  FY11  Tier  1  Warf¬ 
ighter  Outcome.  Enhanced  protection  and  decontamination  of  facilities 
can  be  gained  by  increased  understanding  of  secondary  aerosolization  of 
bioaerosols. 

Thus,  it  is  important  to  continue  to  expand  knowledge  about  particle  re- 
aerosolization,  to  better  understand  bioaerosol  transport  through  build¬ 
ings  and  the  implications  of  secondary  aerosols  in  decontamination  of  af¬ 
fected  facilities.  The  use  of  molecular  models  described  within  this  report 
is  a  promising  method  to  expand  this  area  of  knowledge  and  increase  pro¬ 
tection  and  mitigation  capabilities  for  a  biological  attack. 
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Abbreviations 

Term 

Description 

3-D 

three-dimensional 

ARL 

Army  Research  Laboratory 

BD 

Brownian  dynamics 

C2 

Command  and  Control  Centers 

CB 

chemical-biological 

CC 

computational  chemistry 

CFD 

computational  fluid  dynamics 

CFU 

colony  forming  unit 

CHARMM 

Chemistry  at  FIARvard  Macromolecular  Mechanics;  a  molecular  simulation  program  with 
broad  application  to  many-particle  systems 

CLDI 

concrete-lined  ductile  iron 

COA 

courses  of  action 

COCOM 

Combatant  Command 

CPCM 

polarizable  continuum  model 

cPVC 

chlorinated  polyvinyl  chloride 

DMMP 

dimethyl-methylphosphonate 

DNA 

deoxyribonucleic  acide 

DoD 

Department  of  Defense 

DPD 

diethyl-p-phenylenediamine 

ECBC 

Edgewood  Chemical  and  Biological  Center 

EPANet 

A  public-domain  hydraulic  analysis  package  for  water  supply  networks;  developed  and 
maintained  by  the  USEPA 

ERDC 

Engineer  Research  and  Development  Center 

G-SMD 

grid-steered  molecular  dynamics 

HPC 

heterotrophic  plate  count 

LC/MS 

liquid  chromatography/mass  spectrometry 

MD 

molecular  dynamics 

MOU 

memorandum  of  understanding 

NAMD 

Nanoscale  Molecular  Dynamics 

OP 

organophosphate 

PCM 

Polarizable  Continuum  Model 

PMF 

potential  of  mean  force 

PS 

pipe  segment 

P-S 

phosphorus-sulfur 

PVC 

polyvinyl  chloride 

RNA 

ribonucleic  acid 

SD 

standard  deviation 

Si02 

silica 

TSA 

tryptic  soy  agar 

TST 

Transition  State  Theory 

UIUC 

University  of  Illinois  at  Urbana-Champaign 
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Term 

USEPA 

WISTL 


Description 

US  Environmental  Protection  Agency 
Water  Integrated  Security  Test  Laboratory 
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